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1.0  INTRODUCTION 


The  RACE  (Ring  Acceleration  Experiment)  project  at  Lawrence  Livermore  National  Laboratory 
has  demonstrated  that  plasma  toroids  can  be  accelerated  as  armatures  in  a  coaxial  rail  gun 
configuration  (Ref.  1).  The  magnetic  fields  of  these  toroids  are  similar  to  those  witbm 
spheromaks  (Ref.  2)  except  these  toroids  remain  off  the  axis  of  symmetry  between  a  pair  of 
coaxial  conducting  cylinders  that  serve  as  the  electrodes  for  the  rail  gun.  The  MARAUDER 
(Magnetically  Accelerated  Rings  to  Achieve  Ultra-high  Density  Energy  and  Radiation)  project  at 
the  Weapons  Laboratory  High  Energy  Plasma  Division  (AWX)  will  attempt  to  carry  the  concept 
to  the  megajoule  energy  level.  There  are  many  possible  applications  for  high  energy  toroids. 
They  may  be  used  to  produce  intense  X-ray  radiation  by  interrupting  their  motion  with  a  metal 
plate  or  focusing  cone.  They  may  provide  reliable  armatures  for  plasma  flow  switches  that  can 
be  used  in  pulsed-power  devices  for  current  shaping.  There  are  also  some  exotic  proposals  for 
accelerated  toroids  including  possibilities  for  fusion  (Ref.  3). 

The  full  MARAUDER  experiment  will  be  a  three-stage  process.  The  first  stage  is  the  formation 
of  the  torus  and  is  performed  with  a  magnetized  plasma  gun.  This  is  similar  to  the  spheromak 
formation  that  is  described  by  Turner  (Ref.  2).  The  second  stage  compresses  the  torus  to  a 
smaller  diameter,  increasing  the  plasma  density  and  magnetic  induction.  The  final  stage  is  axial 
acceleration  of  the  toroid.  The  three  stages  are  represented  schematically  in  Figure  1. 

This  report  describes  an  analytic  model  of  the  compression  stage.  The  model  is  presently  being 
used  to  supplement  time-dependent  two-dimensional  magneto-hydrodynamic  (MHD) 
simulations  that  are  performed  on  computers.  The  information  from  the  two  approaches  will  be 
used  in  the  design  of  the  compression  cone  and  acceleration  circuit.  The  analytic  model  gives 
the  magnetic  field  and  associated  magnetic  energy  for  a  torus  before  and  after  compression.  The 
increase  in  torus  energy  plus  the  field  energy  stored  behind  the  torus  is  an  estimate  of  what  the 
acceleration  circuit  must  provide.  Section  2.0  describes  the  solutions  for  the  magnetic  field  and 
integrated  energy.  Section  3.0  presents  the  method  used  for  finding  the  eigenvalue  for  the 
MARAUDER  geometry,  and  Section  4.0  derives  a  model  from  which  comes  the  energy  estimate. 
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Formation  Compression  Acceleration 

Figure  1.  Schematic  of  the  poloidal  magnetic  flux  during  the  three  stages  of  the  MARAUDER 
experiment. 


2.0  ANALYTIC  SOLUTION 


In  low-p  plasmas  where  hydrostatic  pressure  is  small  compared  to  magnetic  field  pressure,  a 
reasonable  assumption  is  that  the  only  forces  of  interest  are  those  from  the  magnetic  fields. 

These  plasmas  are  "force-free"  when  the  plasma  current  density  is  nonexistent  or  everywhere 
parallel  to  the  magnetic  induction.  The  latter  condition  may  be  expressed  as 

Vx  B=X(  x)  B  (1) 

where  B  is  the  magnetic  induction  vector.  The  only  restriction  on  this  equilibrium  configuration 
is  that  B  •  VA.(  x )  =  0  to  enforce  V  •  B  =  0.  Woltjer  (Ref.  3)  has  shown  that  the  integral  of  the 
dot  product  of  B  and  A,  the  magnetic  vector  potential,  over  the  volume  of  a  closed  system  is 
conserved  when  the  system  obeys  the  ideal  MHD  approximation.  This  integral  is  now  known  as 
the  magnetic  helicity,  and  Woltjer  also  proved  that  minimizing  the  magnetic  energy  of  such  a 
system  keeping  its  helicity  constant  leads  to  Equation  1  with  constant  A..  A  minimum  energy 
configuration  is  stable,  and  the  fields  within  spheromaks  and  the  toroids  of  RACE  and 
MARAUDER  have  such  configurations.  It  is  also  useful  to  note  that  A  satisfies 

_  _  _  _  3 

Vx  tfk  +  V<j>,  where  the  gauge  determines  <j>,  so  the  magnetic  energy,  (l/2p-o>J  B-  B  d  x,  is 
A/2  times  the  helicity  (Ref.  4). 

Finding  the  solution  of  Equation  1  with  constant  A.  subject  to  appropriate  boundary  conditions  is 
not  new  (Refs.  4  and  5,  for  example).  The  method  is  presented  here  for  completeness.  The 
energy  density  of  the  resulting  fields  is  then  integrated  for  the  compression  analysis.  The 
conventions  of  Finn,  Manheimer  and  Ott  (Ref.  4)  are  used  except  cylindrical  symmetry  is 
assumed,  so  the  scalar  variations  with  respect  to  the  azimuthal  coordinate  are  zero.  This 
assumption  implies  but  does  not  prove  that  a  symmetric  configuration  has  a  smaller  A.,  hence 
lower  energy,  than  any  other  configuration. 

The  solution  of  Equation  1  with  constant  A.  may  be  obtained  from  the  solution  of  the  scalar 
Helmholtz  equation, 

V2\jf  +  A.\y  =  0  (2) 

The  magnetic  induction  vector  is  then  determined  by 
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B  =  z  x  Vy  +  (1A)V  x  (  z  x  Vy ) 


(3) 


where  z  is  the  axial  unit  vector  in  cylindrical  coordinates.  Morse  and  Feshbach  (Ref.  6)  cover 
this  technique  in  detail  for  generalized  vector  fields.  The  domain  of  this  model  is  an  annular 
volume  that  forms  a  rectangle  when  sliced  by  the  r-z  plane--see  Figure  2.  The  rectangle  has  two 


axis 

Figure  2.  The  rectangular  annulus  that  forms  the  domain  of  the  analytic  model  and  its 
intersection  with  the  r-z  plane. 

sides  parallel  to  z.  The  boundary  conditions  that  make  B  tangent  to  the  bounding  surfaces  are 
2 

d  y 

- =  0  (4) 

5rdz 

at  r  =  r  and  r  =  r  ,  and 
1  2 

2 

9y  2 

—  +4  =  0  (5) 

dz 

at  z=0  and  z=L,  where  r  is  the  radial  coordinate,  r(  and  r^  are  the  inner  and  outer  radii  of  the 


formation  section  (see  Figure  1),  and  L  is  the  axial  length.  Equation  2  may  be  solved  by 
separation  of  variables,  setting  y(r,z)  =  R(r)Z(z),  and  the  differential  equation  becomes 


i  df  aR^i  2  1  a2z 

rR  dr  9r  J  Z  d  z 


2 

The  two  sides  of  the  equation  must  be  equal  to  a  constant  a  .  The  resulting  ordinary  differential 
equation  for  Z  has  a  solution  that  is  the  sum  of  a  sine  function  and  a  cosine  function. 


Z  =  C  sin(az)  +  C  cos(ocz) 

s  c 

The  boundary  conditions  then  require 

C  (X  -  a  )sin(az)  +  C  (X  -  a  )cos(otz)  =  0 

s  c 


at  z  =  0  ,  z  =  L  .  To  avoid  restricting  the  rest  of  the  solution,  the  condition  at  z  =  0  requires 
C  =  0  .  A  series  of  eigenvalues,  a  =  nrc/L  with  n  =  1,  2, ...,  for  a  series  of  sine  terms,  will  then 

c  n 

satisfy  the  boundary  conditions. 


The  solution  of  the  differential  equation  for  R  is  a  sum  of  linearly  independent  Bessel  functions 
of  order  zero.  Before  applying  the  boundary  conditions  at  the  radial  boundaries,  the  solution  for 
y  is 


V(r,z)=  !C[J(r 

n=l  n  0 


a  )  +  A  Y 


(rf^-a2 


:  )]sin(a  z) 

n  n 


If  r  =  0  to  create  a  spheromak,  then  the  boundary  condition  at  r;  forces  A  =  0  for  all  n.  This  is 
not  appropriate  for  coaxial  configurations.  Substituting  Equation  6  into  Equation  4  and 
evaluating  at  the  two  radial  boundaries  gives  two  series  that  are  each  zero.  The  differential 
equation  for  the  minimum  energy  configuration  has  a  constant  X,  and  this  can  be  satisfied  when 
C  is  nonzero  for  only  one  n.  Because  the  energy  is  the  helicity,  which  is  fixed,  multiplied  by 

n 

X/2,  the  configuration  with  the  smallest  X  has  the  lowest  energy.  The  arguments  of  the  Bessel 
functions  must  be  real  for  the  cylindrically  symmetric  case,  as  the  modified  Bessel  functions  that 
would  arise  from  imaginary  arguments  are  monotonic  and  would  not  satisfy  Equation  4  at  both 
radii.  Therefore,  X  >  a  .  Solving  the  n=l  terms  of  the  two  series  with  a  =  n/L, 
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(7a) 


J  (r 
r  1 


yf; 


X  -  a2)  +  AY  (r 

'  1  i 


i, 


X  -  a  )  =  0 


and 


J  (r 

r  2 


Vx2  -  a)  +  AY  (r  ^  -  a2)  =  0 


1  2 


(7b) 


for  the  smallest  radial  factor,  sX  -  a  ,  and  corresponding  A  gives  the  solution  with  the  minimum 
energy.  Determining  the  radial  factor  and  A  is  the  subject  of  Section  3.0. 


Substituting  Equation  6  into  Equation  3  gives  the  components  of  B.  The  constant  is  changed  for 
convenience  to  K  =  C(  Jx  -  a  ,  and  the  components  are 


B  =  (Ka/X)cos(az)[J  ( 

r  1 


-  a)  +  AYt( 


-a2)] 


(8a) 


n  2 

B0  =  -  Ksin(az)[Ji(rvX  -  a  )  +  AY}( 


2 

-a )] 


(8b) 


B  =  -  (Ka 

Z 


J, 


X  -  a  A-)sin(az)(J 


(rJx2  -  a.2)  +  AY  (rlx  -  a 

'  o 


)J 


(8c) 


2 

The  energy  of  this  configuration  is  the  volume  integral  of  the  magnetic  energy  density,  B  /2p.o. 
This  yields 


E 

B 


2  2 .  2  2 

K  Jt  OC  +  X  22  22 

- \ - 2-7 - r[(x  /2)F  (x)  -  xF(x)G(x)  +  (x  /2)G  (x)] 

2\io  {  ctX  (X  -  a) 


■  (1/oX2)[(x2/2)F2(x)  +  (x2/2)G2(x)]  | 


(9) 


where 


F(x)  =  JQ(x)  +  A  Yq(x) 


and 
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G(x)=Ji(x)  +  AY(x) 

When  A  and  X  satisfy  the  boundary  conditions,  G(x)  is  zero  at  the  limits  of  integration,  and  the 
energy  in  the  toroidal  field  is  half  of  the  total.  The  flux  of  toroidal  field  through  the  r-z  plane  is 
also  necessary  for  the  compression  analysis.  The  area  integral  of  Equation  8b  gives 


<J>  =  K — pr - 2  F(x) 

ayX  -  a 


x  =  r  vX  -  a 


GO) 


a 
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3.0  RADIAL  BOUNDARY  CONDITIONS 

The  two  unknowns  in  Equations  7a  and  7b  are  determined  with  a  numerical  root-solver  using 
SLATEC  library  routines  for  the  Bessel  functions.  First,  a  cost  function  is  defined  as  the  sum  of 
the  squares  of  the  residuals  from  Equations  7a  and  7b  when  guesses  are  substituted  for  A  and  X. 
When  the  boundary  conditions  are  satisfied,  the  cost  function  is  a  local  minimum;  in  fact  it  will 
be  zero.  Some  knowledge  of  the  Bessel  functions  and  the  desired  results  helps.  The  Bessel 
functions  are  smooth,  and  the  smallest  radial  factor  is  desired. 

The  program  BESROOT  (see  the  Appendix  for  a  listing)  provides  a  semi-automated  approach  to 
finding  the  roots  in  the  two-dimensional  space  of  A  and  X.  The  user  selects  a  set  X’s,  and  the 
program  finds  the  corresponding  set  of  minimum  costs  between  A  =  -10  and  A  =  +10  using  a 
Newton’s  method  in  the  A  dimension.  A  larger  range  of  A  is  possible,  but  this  would  imply  that 
the  Bessel  function  of  the  second  kind  by  itself  would  almost  fit  the  geometry.  Also,  the  search 
in  X  could  have  been  automated,  but  this  would  add  much  more  complexity,  and  the  manual 
search  is  not  time  consuming  for  one  geometry. 

The  MARAUDER  formation  stage  dictates  r  =  0.4477  m  and  r  =0.6255  m,  the  inner  and  outer 

*  +[2  2 

radii,  respectively,  of  the  chamber.  The  smallest  radial  factor,  \X  -  a  ,  that  satisfies  the 
boundary  conditions  is  17.744  m  ,  and  the  corresponding  A  is  1.3191.  The  axial  length,  L,  of 
the  toroid  is  somewhat  arbitrary.  The  axial  confinement  of  the  plasma  is  a  dynamic  process  in 
the  experiment  and  is  beyond  .ne  scope  of  this  analysis.  Good  confinement  would  produce  a 
toroid  whose  L  is  about  the  same  as  the  radial  gap  between  the  conducting  walls.  Therefore,  for 
one  example,  L  is  set  to  0.2  m,  and  X  becomes  23.698  m  .  Figure  3  shows  the  resulting  field 
distribution  serving  a  second  role  as  initial  conditions  for  a  time-dependent  computer  simulation. 
This  simulation  has  K  =  1,  which  gives  4.388  kJ  of  magnetic  field  energy. 
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L. 


(a)  Poloidal  magnetic  field  vectors. 


(b)  Toroidal  magnetic  field  contours. 

Figure  3.  The  magnetic  fields  of  the  analytic  model  as  initial  conditions  for  a  computer 
simulation. 
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4.0  COMPRESSION  ANALYSIS 


Several  assumptions  and  approximations  are  necessary  to  gain  useful  information  about  the 
compression  stage  from  the  expressions  that  have  been  derived.  The  analysis  provides  the 
energy  of  the  compressed  state,  and  the  assumptions  deal  with  what  the  compressed  state  is 
without  solving  time-dependent  equations.  The  preceding  discussion  has  dealt  with  "force-free" 
magnetic  fields  within  a  plasma-filled  rectangular  annulus  and  has  not  addressed  what  exists 
outside  of  the  annulus.  If  it  were  surrounded  by  a  perfectly  conducting  surface,  eddy  currents 
would  run  along  the  surface.  Material  strength  of  the  surface  would  be  required  to  oppose  the 
pressure  of  the  poloidal  field,  which  does  not  penetrate  the  surface.  If  there  were  no  conducting 
surface,  the  configuration  could  remain  in  equilibrium  only  if  vacuum  fields  matched  the 
poloidal  fields  at  the  surface  of  the  annulus.  The  experimental  toroids  of  MARAUDER 
represent  neither  extreme,  but  they  have  elements  of  each.  The  radial  boundaries  of  the  chamber 
are  only  18  cm  apart  in  the  formation  section,  and  the  total  chamber  length  can  be  several  meters. 
The  concept  is  to  form  a  confined  plasma  ring,  and  accelerate  it  through  the  chamber.  The  goal 
is  not  to  fill  the  entire  chamber  with  plasma.  Vacuum  fields  and  inertia  help  the  axial 
confinement.  The  whole  process  is  dynamic,  and  equilibrium  is  never  established. 


The  compression  process  is  approximated  by  forcing  the  torus  to  remain  in  a  rectangular  annulus 
whose  dimensions  collapse  linearly  towards  a  point  on  the  axis.  Different  acceleration  tube 
geometries  may  be  explored  by  changing  the  final  size  of  the  theoretical  annulus.  The  inner 
radius,  r ,  of  the  annulus  is  a  convenient  parameter,  and  the  outer  radius,  r>,  and  axial  length,  1, 
are  proportional  to  it: 


r  (r. )  =  xr  ,  and  l(r.)  =  qr. 

oil  i  i 


The  constant  q  is  determined  by  the  formation  chamber,  and  X  comes  from  the  arbitrary  choice 
for  L.  For  the  configuration  to  remain  unchanged  with  various  r.,  the  arguments  of  the 
trigonometric  functions  and  Bessel  functions  in  Equations  8a  through  8c  must  remain 
unchanged.  Therefore, 


a(r )  =  rc/qr  ,  and  X(r.)  =  A/r 

it  it 


where  A  =  r^Xfr ).  It  is  also  convenient  to  evaluate  the  definite  integral  in  the  energy  expression 
that  satisfies  the  boundary  conditions,  because  it  will  also  remain  unchanged  for  various  r.: 
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X  2 

n=— tF(x)] 
2 


Its  value  is  2.763  for  MARAUDER. 


If  the  plasma  obeys  ideal  MHD,  where  the  plasma  acts  as  a  perfectly  conducting  fluid,  then  the 

magnetic  flux  through  a  surface  element  of  the  plasma  remains  constant.  Thus,  the  flux 

2 

expression.  Equation  10,  implies  K(r )  =  K(r  )(r /r  ) .  The  magnetic  energy  in  the  torus  is  now  a 
function  of  r  only. 


2  2  2 
H0( A  -  it  h\  )r. 


(ID 


The  magnetic  energy  of  the  idealized  torus  is  inversely  proportional  to  radius.  A  compression 
cone  that  has  been  proposed  for  MARAUDER  reduces  the  radius  by  a  factor  of  three,  so  the 
torus  illustrated  in  Figure  3  would  emerge  from  the  cone  with  approximately  13  kJ  of  energy. 


The  acceleration  circuit  provides  azimuthal  magnetic  field  behind  the  torus  to  compress  and 
accelerate  it.  The  current  in  the  acceleration  circuit  peaks  when  the  torus  reaches  the  end  of  the 
compression  cone.  The  circuit  is  crowbarred  at  this  point  to  efficiently  accelerate  the  torus  as  the 
stored  field  expands  adiabatically.  An  estimate  of  the  energy  the  acceleration  circuit  must 
provide,  E^,  for  a  given  torus  is  the  work  done  on  the  torus  during  compression  plus  the  field 
energy  stored  behind  the  torus, 

E  =  AE  +  (1/2)(L  +  L  )I2  (12) 

c  B  in  ex 

The  internal  energy,  kinetic  energy  and  radiative  losses  are  relatively  small  prior  to  acceleration. 
The  inductance  per  unit  length  between  concentric  cones  is  the  same  as  that  between  coaxial 
cylinders  with  the  same  ratio  of  radii,  %.  Therefore  the  inductance  internal  to  the  chamber  is 


=  (p.oh/2rc)ln(x) 

where  h  is  the  axial  distance  between  the  bottom  of  the  formation  section  and  the  end  of  the 
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compression  cone.  The  inductance  external  to  the  chamber,  L^,  is  determined  experimentally  by 
shorted  discharges. 


One  final  approximation  is  necessary  to  relate  the  peak  current,  I,  to  the  compressed  torus.  If  the 
torus  at  the  end  of  the  cone  is  in  equilibrium  with  the  field  from  the  circuit,  then  the  force  from 
one  balances  that  of  the  other.  These  forces  do  not  have  the  same  distribution,  and  the  circuit 
field  will  distort  the  upstream  side  of  the  torus.  However,  if  an  imaginary  conducting  washer 
separates  the  torus  from  the  circuit  field,  then  the  area  integrals  of  magnetic  pressure  can  be 
matched  to  provide  an  expression  for  I.  The  toroid  has  only  radial  field  at  the  axial  boundaries 
of  the  rectangular  annulus,  and  using  Equation  8a  its  force  on  the  washer  is 


[K^)]2 


3  4 

7t  Or 

l 


2  2  2  2  2  2 
H0T1  A  (A  -  71  /n  )r. 


and  the  force  from  the  circuit  is  (n  I  /47t)ln(x),  so  the  current  that  balances  a  torus  is 


2  2 
2  7t  r 


I=K(r) 


1  ^oTi 


Q 


,1/2 


"2  T 


ln(x)(A  -  7t  A]  ) 


The  complete  expression  for  the  required  circuit  energy  is 


E  = 


C 


[K^Wr,3 

2  2  T~ 

n0(A  -  7t  /n  ) 


A 

j 


ja  h 
271 


y 


[K(ri)]247t4Qri2 
T~ 1  2  2  2  2 

M-0  A  (A  -  7t  (t\  )ln(x) 


J 


(13) 


This  is  proportional  to  the  square  of  the  initial  field,  through  the  K(r  )  factor,  and  it  is  a  quadratic 
of  the  compression  ratio,  r /r.,  where  the  highest  order  term  is  the  stored  field  energy  behind  the 
torus.  Finally  the  length  of  the  cone  is  a  linear  multiplier  in  the  internal  inductance  term.  With 
h  =  0.7  m,  and  L  =50  nH--one  possible  combination  for  MARAUDER,  the  peak  circuit  current 

ex 

for  the  toroid  in  Figure  3  is  1.610  MA,  and  the  circuit  must  provide  a  total  of  134  kJ. 


12 


5.0  DISCUSSION 


This  analysis  for  the  compression  of  a  torus  is  not  a  detailed  model.  The  issue  of  axial 
confinement  is  neglected  as  the  imaginary  conducting  surface  in  the  model  always  encloses  the 
plasma  and  magnetic  fields.  This  is  a  gross  approximation  when  one  considers  that  the  same 
configuration  placed  in  the  actual  chamber  (without  surfaces  on  the  axial  boundaries)  does  not 
satisfy  Ampere’s  Law.  Therefore,  the  real  torus  must  have  vacuum  fields  surrounding  it. 
However,  the  energy  expression  from  the  model  does  provide  useful  information.  It  gives  a 
rough  idea  of  what  circuit  parameters  are  legitimate  for  a  given  torus  and  cone  geometry.  The 
analysis  is  not  complete,  but  it  does  help  limit  the  parameter  space  of  the  computational  study, 
saving  time  and  computer  resources. 
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APPENDIX 


program  besroot 

c - this  program  helps  find  the  solution  to  the  equation 

c  Jl(x)  +  a*Yl(x)  =  0,  where  x  =  r*sqrt(mu**2-(pi/L)**2), 
c  at  two  values  of  r. 

real  sinev,  besev,  totev,  evinc,  rl,  r2,  a,  amx,  amn 
integer  ita,  ite 

common  /radii/  rl,  r2 

open(  unit  =  1,  status  =  ’new’,  file  =  ’besout’) 
call  link(’unit2=tty//’) 

c - for  L  =  20  cm 

sinev  =  15.708 

c - formation  region 

rl  =  0.4477 
r2  =  0.6255 
amaxi  =  10. 
amini  =  -10. 
ita  =  40 


100  write(*,*)  ’How  many  iterations  of  the  eigenvalue,  ’ 
write(*,*)  ’what  starting  point,  and  increment?' 
read(2,*)  ite,  totev,  evinc 

do  200  il  =  l,ite 
amx  =  amaxi 
amn  =  amini 
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besev  =  sqrt(  totev*  *2  -  sinev**2  ) 
dcostmax  =  dcost(besev,  amx) 
dcostmin  =  dcost(besev,  amn) 
if  ( (dcostmax  *  dcostmin)  .It.  0. )  then 

c - the  pair  stradle  0. 

do  300  i2  =  l,ita 
an  =  (amx  +  amn)  /  2. 
dcostn  =  dcost(besev,  an) 

* 

if  ( (dcostmax  *  dcostn)  .It.  0. )  then 
amn  =  an 

dcostmin  =  dcost(besev,  amn) 
else 

amx  =  an 

dcostmax  =  dcost(besev,  amx) 
endif 

300  continue 

write(*,*)  totev,  amn,  valu(besev,amn,rl),  valu(besev,amn,r2) 
write(l,*)  totev,  amn,  valu(besev,amn,rl),  valu(besev,amn,r2) 
else 

write(*,*)  totev,  ’No  root  between’,  amn,  amx 
write(l,*)  totev,  ’No  root  between’,  amn,  amx 
endif 

totev  =  totev  +  evinc 
200  continue 

write(*,*)  ’Type  l  for  another  round.’ 
read(2,*)  irep 
if  (irep  .eq.  1)  goto  100 
close(unit  =1) 

stop 

end 
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function  dcost(ev,  var) 


c - derivative  of  the  cost  function  with  respect  to  var 

c  where  the  cost  is  the  sum  of  the  squares  of  the  residual 
c  of  the  function  in  the  first  comment 

common  /radii/  r  1 ,  r2 

besjll  =  besjl(  ev*rl ) 
besjl2  =  besjl(  ev*r2 ) 
besyll  =  besyl(  ev*rl ) 
besyl2  =  besyl(  ev*r2 ) 

dcost  =  2.  *  ( besjll  +  var  *  besyll )  *  besyll  + 

%  2.  *  (  besjl2  +  var  *  besyl2 )  *  besyl2 

return 

end 

function  valu(ev,  var,  r) 

c - evaluate  the  function  in  the  first  comment 

arg  =  ev  *  r 

valu  =  besjl(arg)  +  var  *  besyl(arg) 

return 

end 
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